*
* $Id$
*
* $Log: twob.F,v $
* Revision 1.2  2002/07/10 09:45:00  morsch
* Gheisha corrections suggested by Gary Bower (FNAL).
*
* Revision 1.1.1.1  1995/10/24 10:20:59  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
*-- Author :
      SUBROUTINE TWOB(IPPP,NFL,AVERN)
C
C *** GENERATION OF MOMENTA FOR ELAST. AND QUASI ELAST. 2 BODY REACT. ***
C *** NVE 04-MAY-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT 15-SEP-1987
C
C THE SIMPLE FORMULA DS/D|T| = S0* EXP(-B*|T|) IS USED.
C THE B VALUES ARE PARAMETRIZATIONS FROM EXPERIMENTAL DATA .
C NOT AVAILABLE VALUES ARE TAKEN FROM THOSE OF SIMILAR REACTIONS
C
#include "geant321/s_defcom.inc"
      DIMENSION RNDM(3)
C
C
      REAL ZNOL,ATNOL
C
C     DATA CB/3./
      DATA CB/0.01/
C
C --- STATEMENT FUNCTIONS ---
      BPP(X)=4.225+1.795*LOG(X)
C
C**
C**  FOR DIFFRACTION SCATTERING ON HEAVY NUCLEI USE BETTER ROUTINE
C**  "COSCAT"
C
      ZNOL = ZNO2
      IF (NFL .EQ. 1)ZNOL = ZNOL - 1
      ATNOL = ATNO2 - 1
      TARMAS=RMASS(14)
      IF (NFL .EQ. 2) TARMAS=RMASS(16)
      ENP(8)=RMASS(IPPP)**2+TARMAS**2+2.0*TARMAS*ENP(6)
      ENP(9)=SQRT(ENP(8))
      EK=ENP(5)
      EN=ENP(6)
      P=ENP(7)
      S=ENP(8)
      RS=ENP(9)
      CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
      IF(ATNO2.LT.1.5) GOTO 500
      IPA1=ABS(IPA(1))
      IPA2=ABS(IPA(2))
      RMC=RMASS(IPA1)
      RMD=RMASS(IPA2)
      RCHC=RCHARG(IPA1)
      RCHD=RCHARG(IPA2)
      IF(ABS(RMC-AMAS).GT.0.001) GOTO 500
      RMNVE=RMASS(14)
      IF (NFL .EQ. 2) RMNVE=RMASS(16)
      IF(ABS(RMD-RMNVE).GT.0.001) GOTO 500
      IF(ABS(RCHC-NCH).GT.0.5) GOTO 500
      IF(NFL.EQ.1.AND.RCHD.LT.0.5) GOTO 500
      IF(NFL.EQ.2.AND.ABS(RCHD).GT.0.5) GOTO 500
      IF(ENP(1).GT.0.0001.OR.ENP(3).GT.0.0001) GOTO 500
      CALL COSCAT
      GO TO 9999
C**
C**  SET EFFECTIVE 4-MOMENTUM OF INITIAL PARTICLE
C**
  500 PV( 1,MXGKPV-1)=P*PX
      PV( 2,MXGKPV-1)=P*PY
      PV( 3,MXGKPV-1)=P*PZ
      PV( 4,MXGKPV-1)=EN
      PV( 5,MXGKPV-1)=AMAS
      PV( 6,MXGKPV-1)=NCH
      PV( 7,MXGKPV-1)=TOF
      PV( 8,MXGKPV-1)=IPART
      PV( 9,MXGKPV-1)=0.
      PV(10,MXGKPV-1)=USERW
      IER(47)=IER(47)+1
      IF(NPRT(4))
     $  WRITE(NEWBCD,4001) (PV(J,MXGKPV-1),J=1,10),IPA(1),IPA(2)
      DO 2 J=1,6
    2 PV(J,1)=PV(J,MXGKPV-1)
      PV(7,1)=1.
      IF(PV(5,1).LT.0.) PV(7,1)=-1.
      PV(5,1)=ABS(PV(5,1))
      NT=1
C**
C** TWO-BODY SCATTERING POSSIBLE?? IF NOT, CONTINUE WITH ORIGINAL
C** PARTICLE, BUT SPEND THE NUCLEAR EVAPORATION ENERGY
C**
      IF(P.LT.0.1) GOTO 200
      IF(RS.LT.0.01) GOTO 200
C**
C** CALCULATE SLOPE B FOR ELASTIC SCATTERING ON PROTON/NEUTRON
C**
      B=BPP(P)
      IF(B.LT.CB) B=CB
      IF(ABS(IPA(2)).GT.13) GOTO 9
      IPA(2)=14
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.5) IPA(2)=16
C**
C** SET MASSES AND MOMENTA FOR FINAL STATE PARTICLES
C**
    9 RMC=RMASS(ABS(IPA(1)))
      RMD=RMASS(ABS(IPA(2)))
      PV(6,1)=RCHARG(ABS(IPA(1)))
      PV(6,2)=RCHARG(ABS(IPA(2)))
      PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD
      IF(NPRT(4)) WRITE(NEWBCD,4002) RMC,RMD,PV(6,1),PV(6,2),RS,S,PF
      IF(PF.LT.0.001) GO TO 9999
      PF=SQRT(PF)/(2.*RS)
C**
C** SET BEAM AND TARGET IN CMS
C**
      PV(1,3)=0.
      PV(2,3)=0.
      PV(3,3)=P
      PV(5,3)=ABS(AMAS)
      PV(4,3)=SQRT(P*P+AMAS*AMAS)
      PV(1,4)=0.
      PV(2,4)=0.
      PV(3,4)=0.
      RMNVE=RMASS(14)
      IF (NFL .EQ. 2) RMNVE=RMASS(16)
      PV(4,4)=RMNVE
      PV(5,4)=RMNVE
C**
C** TRANSFORM INTO CMS.
C**
      CALL ADD(3,4,10)
      CALL LOR(3,10,3)
      CALL LOR(4,10,4)
C**
C** SET FINAL STATE MASSES AND ENERGIES IN CMS
C**
      PV(5,1)=ABS(RMC)
      PV(5,2)=ABS(RMD)
      PV(7,1)=1.
      PV(7,2)=1.
      IF(RMC.LT.0.) PV(7,1)=-1.
      IF(RMD.LT.0.) PV(7,2)=-1.
      PV(4,1)=SQRT(PF*PF+PV(5,1)*PV(5,1))
      PV(4,2)=SQRT(PF*PF+PV(5,2)*PV(5,2))
C**
C** SET |T| AND |TMIN|
C**
      CALL GRNDM(RNDM,2)
      CALL LENGTX(3,PIN)
      BTRANG=B*4.*PIN*PF
C**
C** SIMPLY A PROTECTION AGAINST EXPONENT OVERFLOW 1.E20 IS BIG ENOUGH
C**
      EXINDT=-1.
      IF(BTRANG.LT.46) EXINDT=EXINDT+EXP(-BTRANG)
      TDN=LOG(1.+RNDM(1)*EXINDT)/BTRANG
C**
C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI
C**
      CTET=1.+2.*TDN
      IF(ABS(CTET).GT.1.) CTET=SIGN(1.,CTET)
      STET=SQRT((1.-CTET)*(1.+CTET))
      PHI=RNDM(2)*TWPI
C**
C** CALCULATE FINAL STATE MOMENTA IN CMS
C**
      PV(1,1)=PF*STET*SIN(PHI)
      PV(2,1)=PF*STET*COS(PHI)
      PV(3,1)=PF*CTET
      PV(1,2)=-PV(1,1)
      PV(2,2)=-PV(2,1)
      PV(3,2)=-PV(3,1)
C**
C** TRANSFORM INTO LAB
C**
      DO 11 I=1,2
      CALL LOR(I,4,I)
      CALL DEFS1(I,MXGKPV-1,I)
      IF(ATNO2.LT.1.5) GOTO 11
      CALL LENGTX(I,PP)
      IF(PP.LT.0.001) GOTO 11
      EKIN=PV(4,I)-ABS(PV(5,I))
      CALL NORMAL(RAN)
      EKIN=EKIN-CFA*(1.+0.5*RAN)
      IF(EKIN.LT.0.0001) EKIN=0.0001
      PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
      PV(4,I)=EKIN+ABS(PV(5,I))
      PV(1,I)=PV(1,I)*PP1/PP
      PV(2,I)=PV(2,I)*PP1/PP
      PV(3,I)=PV(3,I)*PP1/PP
   11 CONTINUE
      NT=2
C**
C** ADD BLACK TRACK PARTICLES .
C** HERE THE PROCEDURE IS SOMEWHAT DIFFERENT AS IN 'TWOCLU' AND 'GENXPT'
C** THE REASON IS, THAT WE HAVE TO SIMULATE ALSO THE NUCLEAR REACTIONS
C** AT LOW ENERGIES LIKE A(H,P)B, A(H,P P)B, A(H,N)B E.T.C.
C**
  200 IF(ENP(1).LE.0.0001.AND.ENP(3).LE.0.0001) GOTO 40
      IF(ATNOL .LT. 0.9)GOTO 445
      SPALL=0.
      TEX=ENP(1)
      IF(TEX.LT.0.0001) GOTO 445
      BLACK=TEX/0.02
      CALL POISSO(BLACK,NBL)
      IF(NBL.GT.ATNOL) NBL=ATNOL
      IF(ENP(1).GT.0.0001.AND.NBL.LE.0) NBL=1
      IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX
      IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
      IF(NBL.LE.0) GOTO 445
      EKIN=TEX/NBL
      EKIN2=0.
      CALL STEEP(XX)
      DO 441 I=1,NBL
      IF(NT.EQ.MXGKPV-2) GOTO 441
      IF(EKIN2.GT.TEX) GOTO 443
      CALL GRNDM(RNDM,1)
      RAN1=RNDM(1)
      CALL NORMAL(RAN2)
      EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
      IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
      EKIN1=EKIN1*XX
      EKIN2=EKIN2+EKIN1
      IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
      IF(EKIN1.LT.0.) EKIN1=0.0001
      IPA1=16
      PNRAT=1.-ZNOL/ATNOL
      CALL GRNDM(RNDM,3)
      IF(RNDM(1).GT.PNRAT)THEN
      IPA1=14
      ZNOL = ZNOL -1
      ENDIF
      ATNOL = ATNOL -1
      NT=NT+1
      SPALL=SPALL+1.
      COST=-1.+RNDM(2)*2.
      SINT=SQRT(ABS(1.-COST*COST))
      PHI=TWPI*RNDM(3)
      IPA(NT)=-IPA1
      PV(5,NT)=ABS(RMASS(IPA1))
      PV(6,NT)=RCHARG(IPA1)
      PV(7,NT)=2.
      PV(4,NT)=EKIN1+PV(5,NT)
      PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
      PV(1,NT)=PP*SINT*SIN(PHI)
      PV(2,NT)=PP*SINT*COS(PHI)
      PV(3,NT)=PP*COST
  441 CONTINUE
  443 IF(ATNO2.LT.10.) GOTO 445
      IF(EK.GT.2.0) GOTO 445
      II=NT+1
      KK=0
      EKA=EK
      IF(EKA.GT.1.) EKA=EKA*EKA
      IF(EKA.LT.0.1) EKA=0.1
      IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
      IF(IKA.LE.0) GO TO 445
      DO 444 I=1,NT
      II=II-1
      IF(IPA(II).NE.-14) GOTO 444
      IPA(II)=-16
      IPA1  = 16
      PV(5,II)=ABS(RMASS(IPA1))
      PV(6,II)=RCHARG(IPA1)
      KK=KK+1
      IF(KK.GT.IKA) GOTO 445
  444 CONTINUE
  445 TEX=ENP(3)
      IF(TEX.LT.0.0001) GOTO 40
      NBL=IFIX(2.*LOG(ATNO2))
      IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
      IF(NBL.LE.0) GOTO 40
      IF(ATNOL .LT. ZNOL + 1)GOTO 40
      IF(ZNOL .LT. 1.0)GOTO 40
      EKIN=TEX/NBL
      EKIN2=0.
      CALL STEEP(XX)
      IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX
      DO 442 I=1,NBL
      IF(ATNOL .LT. ZNOL + 1)GOTO 40
      IF(ZNOL .LT. 1.0)GOTO 40
      IF(NT.EQ.MXGKPV-2) GOTO 442
      IF(EKIN2.GT.TEX) GOTO 40
      CALL GRNDM(RNDM,1)
      RAN1=RNDM(1)
      CALL NORMAL(RAN2)
      EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
      IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
      EKIN1=EKIN1*XX
      EKIN2=EKIN2+EKIN1
      IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
      IF(EKIN1.LT.0.) EKIN1=0.0001
      CALL GRNDM(RNDM,3)
      COST=-1.+RNDM(1)*2.
      SINT=SQRT(ABS(1.-COST*COST))
      PHI=TWPI*RNDM(2)
      RAN=RNDM(3)
      IPA(NT+1)=-30
      ATNOL = ATNOL -2
      ZNOL = ZNOL -1
      IF(RAN.GT.0.60)THEN
      IF(ATNOL .GT. ZNOL + 0.9)THEN
      ATNOL = ATNOL -1
      IPA(NT+1)=-31
      IF(RAN.GT.0.90)THEN
      IF( (ATNOL .GT. 0.9) .AND. (ZNOL .GT. 0.9) )THEN
      ATNOL = ATNOL -1
      ZNOL = ZNOL -1
      IPA(NT+1)=-32
      ENDIF
      ENDIF
      ENDIF
      ENDIF
      INVE=ABS(IPA(NT+1))
      PV(5,NT+1)=RMASS(INVE)
      SPALL=SPALL+PV(5,NT+1)*1.066
      NT=NT+1
      PV(6,NT)=RCHARG(INVE)
      PV(7,NT)=2.
      PV(4,NT)=PV(5,NT)+EKIN1
      PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
      PV(1,NT)=PP*SINT*SIN(PHI)
      PV(2,NT)=PP*SINT*COS(PHI)
      PV(3,NT)=PP*COST
  442 CONTINUE
C**
C** STORE ON EVENT COMMON
C**
   40 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
      EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
      EKIN2=0.
      CALL TDELAY(TOF1)
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      TOF=TOF-TOF1*LOG(RAN)
      DO 1 I=1,NT
      EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
      IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
      PV(7,I)=TOF
      PV(8,I)=ABS(IPA(I))
      PV(9,I)=0.
    1 PV(10,I)=0.
      IF (NPRT(4)) WRITE(NEWBCD,1003) NT,EKIN,EKIN1,EKIN2
      INTCT=INTCT+1.
      CALL SETCUR(NT)
      NTK=NTK+1
      IF(NT.EQ.1) GO TO 9999
      DO 50 II=2,NT
      I=II-1
      IF(NTOT.LT.NSIZE/12) GOTO 43
      RETURN
   43 CALL SETTRK(I)
   50 CONTINUE
C
 1002 FORMAT(' *TWOB* ',5F10.4,10X,5F10.4/1H ,7X,5F10.4,10X,5F10.4/
     $ ' LAB SYSTEM FINAL STATE FOUR VECTORS')
 1003 FORMAT(' *TWOB* COMPARISON',2X,I5,2X,3F10.4)
 4001 FORMAT(' *TWOB* ',10F10.4,2X,2I3)
 4002 FORMAT(' *TWOB* ',7F10.4)
 3003 FORMAT(' *TWOB* ',I3,' BLACK TRACK PARTICLES PRODUCED',
     $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
 3004 FORMAT(' *TWOB* ',I5,' HEAVY FRAGMENTS PRODUCED',
     $ ' WITH TOTAL ENERGY OF',F8.4,' GEV')
C
 9999 CONTINUE
C
      END
